//
//  main.cpp
//  opencv_project
//
//  Created by 黄琛 on 2018/12/27.
//  Copyright © 2018 黄琛. All rights reserved.
//

#include "opencv2/opencv.hpp"

using namespace cv;
using namespace std;

int main(void)
{
    Mat image = imread("rice.png");
    
    Mat gray, bw;  // 二值化后的图像
    
    // 由于threshold只支持灰度图像，所以先转换
    cvtColor(image, gray, COLOR_BGR2GRAY);
    // 大津算法阈值化
    threshold(gray, bw, 0, 0xff, CV_THRESH_OTSU);
    
    // 形态学处理，去除噪声
    Mat element = getStructuringElement(MORPH_CROSS, Size(3, 3));
    morphologyEx(bw, bw, MORPH_OPEN, element);
//    morphologyEx(bw, bw, MORPH_CLOSE, element);
    
    // 以下是图像分割
    Mat seg = bw.clone();
    vector<vector<Point>>cnts;
    findContours(seg, cnts, RETR_EXTERNAL, CHAIN_APPROX_SIMPLE);
    
    // 以下进行筛选
    int count = 0;
    Rect rect;
    string strCount;
    double area = 0;
    double length = 0;
    for(long i = cnts.size() - 1; i >= 0; i--)
    {
        vector<Point> c = cnts[i];
        double area_temp = contourArea(c);
        double length_temp = arcLength(c, true);
//        if(area_temp < 10)  // 滤除面积小于10的分割结果：可能是噪声
//            continue;
        area += area_temp; // 累加米粒面积
        length += length_temp; // 累加米粒轮廓长度
        count++;  // 统计米粒数量
        cout << "blob" << i << ":" << area_temp <<", " << length_temp << endl; // 输出米粒轮廓面积和长度
        rect = boundingRect(c);
        // 在原始图像上画出包围矩阵，并给每个矩形标号
        rectangle(image, rect, Scalar(0, 0, 0xff), 1);
        
        stringstream ss;
        ss << count;
        ss >> strCount;
        putText(image, strCount, Point(rect.x, rect.y), CV_FONT_HERSHEY_PLAIN, 0.5, Scalar(0, 0xff, 0));
    }
    
    double avg_area, avg_length;
    avg_area = area / count;  // 米粒平均面积
    avg_length = length / count;  // 米粒平均轮廓长度
    cout << "米粒数量:" << count << endl;
    cout << "米粒平均面积:" << avg_area << endl;
    cout << "米粒平均长度:" << avg_length << endl;
    
    //计算方差
    double accum_area = 0;
    double accum_length = 0;
    for_each(begin(cnts), end(cnts), [&](vector<Point> c){
        accum_area += (contourArea(c)- avg_area) * (contourArea(c) - avg_area);
        accum_length += (arcLength(c, true) - avg_length) * (arcLength(c, true) - avg_length);
    });
    
    double var_area = accum_area / count;
    double var_length = accum_length / count;
    cout << "米粒面积的方差:" << var_area << endl;
    cout << "米粒长度的方差:" << var_length << endl;
    
    // 筛选面积和长度落在3sigma范围内的米粒
    int count1 = 0;
    int count2 = 0;
    
    for_each(begin(cnts), end(cnts), [&](vector<Point> c){
        if(contourArea(c) > avg_area - 3*sqrt(var_area) && contourArea(c) < avg_area + 3*sqrt(var_area))
            count1++;
        
        if(arcLength(c, true) > avg_length - 3*sqrt(var_length) && arcLength(c, true) < avg_length + 3*sqrt(var_length))
            count2++;
    });
    
    cout << "面积落在3sigma范围内的米粒数量:" << count1 << endl;
    cout << "长度落在3sigma范围内的米粒数量:" << count2 << endl;
    
    imshow("src", image);
    imshow("result of OTSU", bw);
    
    waitKey(0);
    
    return 0;
}
